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Introduction 


GEODYNII is a conventional batch least-squares differential corrector computer program with 
deterministic models of the physical environment Conventional algorithms have been used to 
process differenced phase and pseudorange data to determine eight-day GPS orbits with several 
meter accuracy (Schenewerk, 1990). However, random physical processes drive the errors 
whose magnitudes prevent improving the GPS orbit accuracy. To improve the orbit accuracy, 
these random processes should be modeled stochastically. The conventional batch least-squares 
algorithm cannot accommodate stochastic models, only a stochastic estimation algorithm is 
suitable, such as a sequential filter/smoother. Also, GEODYNII cannot currently model the 
correlation among data values. Differenced pseudorange, and especially differenced phase, are 
precise data types that can be used to improve the GPS orbit precision (Counselman et al., 1989). 
To overcome these limitations and improve the accuracy of GPS orbits computed using 
GEODYNII, we proposed to develop a sequential stochastic filter/smoother processor by using 
GEODYNII as a type of trajectory preprocessor. Our proposed processor is now completed. It 
contains a correlated double difference range processing capability, first order Gauss Markov 
models for the solar radiation pressure scale coefficient and y-bias acceleration, and a random 
walk model for the tropospheric refraction correction. 

The development approach has been to interface the standard GEODYNII output files 
(measurement partials and variationals) with software modules containing the stochastic 
estimator, the stochastic models, and a double differenced phase range processing routine. Thus, 
no modifications to the original GEODYNII software have been required. A schematic of the 
development is shown in Figure 1. The observational data are edited in the preprocessor and the 
data are passed to GEODYNII as one of its standard data types. A reference orbit is determined 
using GEODYNII as a batch least-squares processor and the GEODYNII measurement partial 
(FTN90) and variational (FTN80, V-matrix) files are generated. These two files along with a 
control statement file and a satellite identification and mass file are passed to the filter/smoother to 
estimate time-varying parameter states at each epoch, improved satellite initial elements, and 
improved estimates of constant parameters. 
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Figure 1 . Flowchart showing the procedure as currently developed 

Background 

The following discussion assumes some familiarity with filtering and smoothing theory as 
developed, for example in Brown (1983) and Gelb (1974). Additional familiarity is assumed with 
the square root information filtering and smoothing algorithm as developed in Bierman (1977) and 
implemented by Swift (1987). Any departures in our implementation from Swift's formulation are 
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explained here in detail. Some topics are expanded here to clarify and to supplement the 
discussion in these earlier references. Particular attention has been focused on showing the 
common foundation of the information and covariance filters. The efficiency of the Householder 
Transformation to compute an equivalent square upper triangular matrix from a larger rectangular 
matrix is discussed. Also, the quantities that define the first-order Gauss Markov and random 
walk models are clearly derived. 

The discrete form of the stochastic state equations are: 

Ax J+1 = ^ Ax j + G(Oj ( 1 ) 

where 

A Xj is the state at time t ; 

is the nonsingular state transition matrix relating the state at t ; 
to the state t ;+1 . 

0)j is the vector of white noise process terms with a nonsingular covariance matrix 
Q, with dimco<dimAx. 

G maps the source white noise process into the state with dimAx. 

The discrete form of the linear measurement model is: 

Zj = AjAXj + v j (2) 

where 

z } is the vector of measurements at time t ; 

A j is the matrix of partial derivatives of the measurement model w.r.t. the state at t ; 
v . is the vector of measurement noise with the covariance . 

The observations are decorrelated and whitened so that = I. This is done without a loss of 
generality. A set of observations with F^, = I can be constructed. The procedure will be given 
later. 

A solution to this problem was first proposed by Kalman in early 1960’s (Kalman, 1960; Kalman 
and Bucy, 1961). A solution for the state and its covariance can be derived by applying Bayes' 
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rule. This derivation can be found in Maybeck (1979). The results are repeated here with a slight 
change in notation. 


tk-py + AjAj]" (3) 

tej = Pj.^'AXj + Ajz ; ] (4) 

The symbol refers to the propagated (predicted) estimate of the state and covariance at t ; . 
The ’ A ' symbol refers to the estimate of the state or covariance after incorporating the 
measurement at t ; . The state Ax, is propagated from t y . to time t ;+1 , using equation (1). The 

covariance P, is propagated to t ;>1 , by 

= (5) 

Notice that the inverse of P J+l „ is required in equations (3) and (4). To avoid the inversion of P, +1 

at each step, a direct propagation of F^., is desired. This can be developed by applying the 
following lemma to equation (5). 

(A + X T Y) _1 = A"' - A _1 X T (l + YA- 1 X T )" 1 YA' 1 (6) 

where 

A = O P<D t 

jjj 

X T =G,Q, 

Y = Gj 

and defining 


M >i = 


(7) 
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( 8 ) 


% = M, +1 -C,GJM, +1 , 


where the gain C 7 , is : 


Cj = M^.G^GjM^Gj + Q; , ]‘ I (9) 

To propagate and update the state equations in terms of Ey and Ey , the state estimates are 
replaced by 


y, = E^Ax, (10) 

9j = f^Ai, (11) 

The state update and propagation equations are: 

Yj = yj + Ajz, (12) 

y , + 1 = [i - C 7 gJ ]a Kt j+l ,tj )% ( 1 3) 


The state estimates can be found at any time by solving equations (10) and/or (11) for Ax, and/or 
A Xj. Equations (8) to (13) are an algorithm to solve the problem defined by equations (1) and 

(2). Here, the inverse covariance is propagated. This algorithm is sometimes called a Bayes' 
filter. The inverse covariance is also called an information matrix leading to the name information 
filter. This algorithm requires computing the inverse of an n x n matrix where n is the number of 
states. The state estimate covariance can be completely uncertain since in the inverse of P 0 the 
elements of the matrix become zero. This algorithm is most efficient when the number of 
measurements, m, is relatively larger than the number of states, n, and when the solutions for the 
state and covariance are needed infrequently. 

The usual Kalman filter can be derived from equations (3) and (4) by applying the matrix lemma: 

[F 1 + A T A]" 1 = P- PA t [APA t ] _1 AP (14) 
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These optimal filters, either the Bayes' or Kalman, exhibit numerical instabilities that cause the 
state estimates to diverge (Bierman and Thornton, 1977). More numerically stable gain matrix 
expressions have been derived for both the covariance and information matrix forms (Maybeck 
ibid.). However, these require a significant number of additional matrix computations and are 
thus not completely satisfactory. A more comprehensive approach was to reformulate the filter 
algorithm in terms of square roots of the covariance or the information matrix. The square root 
filter maintains numerical accuracy to approximately the same number of digits with half the word 
length required by a conventional non-square root algorithm. 

The square root (or more correctly the Cholesky factorization) of an nxn matrix N is defined as : 

N=Sg r (15) 

The square roots are not unique. Any orthogonal transformation (T) of a square root matrix S is 
also a square root of N. The useful properties of the orthogonal transformation can be shown by 

factoring I* 1 into the product of its Cholesky factors 




and thus, F^ 1 becomes: 


where 


Now 


E^RJR^RJVAJA, 




(16) 

(17) 

(18) 

(19) 


The transformation T is an orthogonal transformation. Its columns form an orthonormal basis for 
Rj of n vectors since R, has rank of n. The first n vectors span the range space of R y and the 
vectors n+1 to n+m are orthogonal to this spanning set. Thus, the last m rows of R^ are zero. 
Also, the basis vectors were chosen in a manner that R, is upper triangular. A Householder 
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Transformation T is used to compute R ■ (Bierman, ibid.). This allows the square root matrices of 

dimension ( n+m ) by n to be transformed to an equivalent form of an upper triangular square 
matrix of dimension n. 


Square Root Information Filter and Smoother (SRIF/SRIS) 


In application of the Householder Transformation to an augmented matrix is the essence of the 
Square Root Information Filter (SRIF). The equations are most easily constructed using 
Bierman's 'data equation’ point of view. The problem is treated as a least-squares problem where 
the least-squares functional is to minimized. This is accomplished by applying a Householder 
Transformation assuming the state at tj to be a priori information and augmented with the 

measurements at ^ . Thus, 



i 

»N 


• *> 

1 

A 

1 

nT 


i 

o 

& 


( 20 ) 


where the 'data equations' are defined as 

( 21 ) 
( 22 ) 



Swift has shown the equivalence of equation (20) with the more conventional formulation of 
equations (3) and (4). 

The propagation of the state and the covariance were given in equations (1) and (5). These can 
also be incorporated into the SRIF by defining 'data equations' and applying the Householder 
Transformations. The details of the derivations can be found in Bierman (ibid.). The results are 
repeated here. 


■ KU) 

0 

KU) 


A o+D 

R„0 + D I. O' +»1 

-R 

R 4> _1 

jj 

1 

+ 

tfcsT 


0 

J*' 

+ 

s5" 

t 

1 


( 23 ) 


where the 'data equation' for the noise term Q) is 
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z«,0‘) = Ro.OXj) 


(24) 


Swift has also shown the equivalence of equation (23) to equations (1) and (5). 

Bierman further partitions the propagation equation (1) into stochastic states, dynamic states and 
bias states. Equation (1) can be written as: 

Ap M 0 0 Ap 

Ax = V p V x V y Ax + 

.Ay _, +1 L 0 0 I J[Ay| 

where 

Ap is the correlated process noise states 
Ax states that vary with time by not explicitly influenced by process noise. 

Ay bias (constant) parameters. 

Vp, V I} Vj,, are transition matrix elements. 

The dynamic parameters can be redefined in the form of pseudo epoch state parameters. This 
dynamic model definition allows the variationals and measurement partials from a batch 
differential corrector orbit determination program (e.g., GEODYNII) to be used directly in the 
filter algorithms. The state equations now become 

Ap M 0 0 Ap 

Ax = V p I 0 Ax + 

Ay j+1 0 0 1 Ay j 

where 

V„ = v,( W ,) = v' l - | (f J . 1 ,r 0 )v p (( J .„t y ) 

V; l (t J+l ,T 0 ) is the inverse of the state transition matrix interpolated from the 
GEODYNII V-matrix file (FTN80) 

Vp(^,^) is the transition matrix of the time- varying parameters from tj to t J+l . 
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The transformation in equation (20) is written as a two step transformation that saves storing a 
block of ny x (np+nx) zeroes that would be present if equation (20) was used in its original form. 
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The mapping equation (23) becomes, assuming z a is zero 


(28) 
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The subroutine in Bierman ( ibid., 155-157) neglected the upper triangular elements of R p above 
the diagonal. A subroutine to compute right side of equation (29) including the neglected off- 
diagonal elements of R p is given in Appendix B. 


The smoothing process is a backward filter of the forward filter results. For the orbit 
determination problem, the fixed interval smoother is appropriate. For Kalman filtering, the 
Rauch-Tung-Streibel (RTS) smoothing algorithm is widely used (Brown, ibid.). A general 
formulation for inverse covariance smoothing is given by Maybeck (ibid.). The Square Root 
Information Smoother (SRIS) is given by Bierman (ibid.). The equation for the implementation of 
Bierman 's pseudo epoch formulation is given by Swift (ibid.). Swift also shows the equivalence of 
the SRIS to the RTS smoother. The SRIS equation is 
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The top row of the matrix on the right side is not needed again, but the other terms are combined 
with the smoothing coefficients at tj to smooth back to tj_ v 

From the general expression of the data equation, and the relationship of the covariance to the 
square roots of the information matrix, the solution for the states and their covariances for either 
the filter or smoothing operations are 

Xj=R;‘zj (31) 

p y = R;'k; t (32) 

For filtering, the right sides of equations (27) and (28) are solved for the unknowns and the 
covariances as 


~Ap' 



R px 
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(VLfV 

R _1 z 


Ax 

= 
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y y 
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(33) 
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(35) 


p y = r;‘r; t (36) 

For the smoothing problem, the right side of (30) is used in equations (32)-(36) with the 
smoothed value of z y and R r The smoothed values of z y and R y are the values at the last filter 
step. 

Time-Varying Stochastic Parameter Models 

Many physical processes can be modeled using one of the Gauss Markov filters easily constructed 
by filtering white noise through a simple filter. The processes that are of interest here are the 
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first-order Gauss Markov model and the random walk model. The discrete mathematical 
expression for these models is now derived from their continuous forms. The first-order Gauss 
Markov process describes the physical process where the state at t J+l depends only on the 

previous state at t j . This process can be described by the stochastic differential equation. 

^--pit) = (o (37) 

at x 


where 


T is the correlation time 

t o is the white Gaussian noise with zero mean and covariance Q, 

E(q)j ) = 0, E((Oj ,(o k )=Q8 O' - k) = q con 8 (, j - k) where q con is the 
continuous spectral density. 

The variance of p, a p satisfies the differential equation 


da 2 (t) 

dt 


~a 2 p (t) + q a 


(38) 


If the process is allowed to continue for a time interval several multiples longer than t, the 


E(p 2 (f)) will approach a limit and 


dolifi 

dt 


will approach zero. This is the steady state variance. 


dalit) 2 

By setting — 2 — = 0 and solving equation (38) for steady state a p {t st ) is found 

(JLv 


<( 0 = £ Qa 


(39) 


From state space methods the solution of equation (37) is 


j+i 

pij + 1) = M(j + Xj)pij) + J Mij, A )(oiX)dX (40) 

j 

and covariance of p is 


j + l 

C p = MPjAf + J Mij +\X )Q(X)M T ij +\X )dX 


J 


(41) 
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The matrix M is the state transition matrix which must satisfy the relationships 

M=--M and M(j,j) = I 

T 

The solution for M is 

M=e T 

where At = t J+l - 1 } 

Thus the discrete form of the state update becomes 

ApO' + 1) = M(j + \j)Ap(j) + (Oj 
This is the top row of equation (26). 

Now the solution for the discrete covariance update is derived from 


*,♦1 2 U - Q ) 


QdU = | Qcon e <& = ~ 


2 (X-tj) 


Qcon ^ 


Qcon g 


2H \ 


1-e 


The random walk model is a special case of the first order Gauss Markov where r — > °° . 
discrete state update becomes 


Ap(j + 1) = A p(j) + o)j 
The discrete covariance is found from 


Qdi 


= limq ro - 


l-e 
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2At 
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) 


(42) 

(43) 

(44) 

(45) 
The 

(46) 
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= a lim 

'icon _ 


Ar + ^ constant; * j — ) 

i=l 


= Q At 

'icon * 


( 47 ) 


So, to implement the first-order Gauss Markov process, the correlation time (T) and the 
continuous process noise variance (q con ) must be specified. The matrix M is computed using 
equation (43) and the q dil from equation (45). The random walk model is specified by defining 
the continuous process noise variance ( q con ), here M= I, and q dU is computed using equation 
(47). 


Solar Radiation Pressure Scale Coefficient and Y-bias Acceleration Model 

The orbit related stochastic parameters modeled in OSUORBFS are first-order Gauss Markov 
models for the solar radiation pressure scale coefficient and the y-bias acceleration. The V p 

matrix that maps the effects of the stochastic parameters Ap,, on the epoch state parameters Ax^ 
is derived following the scheme of Swift (ibid.). The V p matrix has dimension nx (number of 

pseudo epoch state parameters) by nd (number of orbit-related stochastic parameters) and has the 
general form of 




K 0 
o \ 


0 


3>' 


0 


•. 0 
0 


(48) 


The ith satellite contribution O'*. = (£. +1> f,), a 6 by 2 matrix, to V p is computed as 




ddtj+i) 


'j* i ; 


dK^) <2K,(f,) 

) jgW 




(49) 
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where 


,tj ) is the state transition matrix interpolated from the GEODYNII output 
V-matrix file (FTN80) variationals. See also equation (26). 

The partial derivatives of the position and velocity at </* i with respect to the solar radiation 
pressure scale coefficient and the y-bias acceleration at tj are approximated using a second-order 
Taylor Series expansion. 


where 


ch(t J+l ) _ <2r(^) 

~ 2 cK^itj) 

= 1 = 1,2 
2 dK Ti (tj) ** 


(50) 


®K, -[®K, 



Bk, = 0 for a random walk. 

Since the angle between the x and y satellite axes is not easily estimated it is not modeled in 
OSUORBFS. A 90 degree angle is assumed. Thus the following differs slightly from Swift. 


' Mtj) 

dr(tj) 

= R. 

a™ shape, 
0 

0 

shape 

MW 



a™ shape 

0 


(51) 


where 

a ™ and a™ are the accelerations along the satellite x and z axis respectively. 

The ROCK4 and ROCK42 models are used to compute a"‘ and a™. 
shape is either 0 or 1 depending if the sun is obstructed by the earth from view of 
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the satellite. 

R a is a matrix transformation from the satellite axis system to the True of 
Reference Date (TORD) inertial Cartesian reference system. 

The ROCK4 and ROCK42 models (Fliegel and Gallini, 1992), the matrix transformation R s , and 
the computation of shape require the sun-earth-satellite positions. The mean of date (MOD) 
positions of the sun and earth were computed using the closed form expressions of Fliegel and 
Harrington (1993). 

Tropospheric Refraction Correction 

The measurement-related stochastic parameter modeled in OSUORBFS is a random walk model 
for the refraction correction (Tralli et. al„ 1988; Herring et. al., 1990). This model is defined by 
equations (46) and (47). 

Double Difference Observable Decorrelation and Whitening 

The full covariance matrix for the double difference range data is constructed. The observations 
are then decorrelated and whitened. The general form of the observation equations as defined in 
equation (2) is 


z = Ax + v (52) 

Here, the observation error v has a zero mean, E(v) = 0, but is correlated, E{ w T ) = P v 
A set of uncorrelated observations with unit covariance can be constructed from the lower 
triangular square root of P v . 


P v = L V L V T (53) 

Here, L v can be computed by a lower Cholesky factorization of P v . The desired independent set 
of observations is 


L v ! z = L^Ax + L v V 


(54) 
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At any particular epoch, the m = (#stations-l)x(#satellites-l) linear independent double 
difference range data types can be formed. These m observations are independent in the sense of 
linear algebra, but are statistically correlated. Each of the m observations has the form 
(GEODYNII measurement type=87) 

[(si fl) - (s2 fl)] - [(si -> t2) - (s2 -4 f2)] (55) 


where 

(sl-» fl) etc., are the satellite-station range observations (e.g., satellite 1 to 
station 1). 

For small regional networks, a single satellite station pair is selected as the base satellite-station 
and the m observations are constructed by differencing the remaining satellite-stations with the 
base pair. For a global network, the distance between stations may prevent using a single base 
pair to construct all observations at that epoch. Thus, no consistent numerical structure exists 
that would permit a symbolic construction of the decorrelated measurement set. P v and L" 1 must 
be computed numerically at each epoch. P v is computed using conventional error propagation 

P v = tr r 2 GG T (56) 


where 

cr r 2 is the standard deviation of the single one-way range measurement 
G is the matrix of partial derivatives of the observation equation with respect 

to the one-way range. This matrix contains elements of -1,0,1 which are the linear 
combination of one-way ranges that define the double difference 

The decorrelated observation set with unit covariance is obtain from equations (53) and (54). 

OSUORBFS 

The program OSUORBFS is designed to filter and smooth the GEODYNII batch solutions. 
OSUORBFS requires from GEODYNII the measurement partials file (FTN90) and the 
variationals V-matrix file (FTN80). The user must supply a user input file (FSN05) and a file of 
satellite identification numbers and masses (SATMAS.TAB). The GEODYNII processing 
proceeds in the usual way with TDF, G2S, and G2E program executions. On the last iteration an 
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output of the setup deck (FTN05) which contains the current parameter estimates is requested 
using PUNCH to output the new setup deck in file FTN07. Then FTN07 is modified to include 
global cards PARFIL and EMATRX to output files FTN90 and FTN80. The maximum iteration 
numbers for the global (outer) and arc (inner) are set to one on the ENDGLB and REFSYS 
statements. 

Now TDF, G2S, and G2E are executed and FTN06, FTN80, and FTN90 files are output. An 
alternative approach avoiding the restart of GEODYNII is to force an additional iteration in the 
first GEODYN execution by increasing the outer/inner iteration maximum counts and decreasing 
the RMS tolerance. This approach would be less cumbersome to implemented. Additional 
operational experience is needed to determine which approach is satisfactory. 

The user must construct the control file, FSN05, for OSUORBFS. This file contains six control 
statements (REFSYS, DECORR, FILSMT, UPDTRJ, CONPRT, SATMAS) to control the 
configuration of the filter/smoother solutions. The six statements are mandatory. These 
statements are explained in Appendix A. 

FSN05 must also contain the parameter labels from FTN06. The measurement partials and 
variationals form GEODYNII are identified by internal parameter labels as described in the 
GEODYNII manual volume 5. For OSUORBFS to recognize these partials, the parameter labels 
as they appear in FTN80 and FTN90 must be specified in FSN05. These labels can be accessed 
by printing the EMATRIX header record in FTN06 during the last iteration of the GEODYN run. 
They must be manually edited and placed in FSN05. 

FSN05 must also contain the parameter types as defined in the following table. 


Type# Description 

1 orbit-related stochastic 

2 measurement-related stochastic 

3 pseudoepoch state 

4 measurement-related constant 

5 orbit-related constant 


Example 

solar radiation pressure (1st Gauss Markov) 
y-bias acceleration (1st Gauss Markov) 
tropospheric refraction correction (random 
walk) 

satellite initial elements 

double difference bias, tropospheric refraction 
correction 

solar radiation pressure coefficient, y-bias 
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For each parameter type the apriori standard deviation must be specified. Additionally, for the 
first-order Gauss Markov model, the continuous process noise standard deviation Qq con ) and the 
correlation (T ) time must be specified. For the random walk model, the continuous process noise 
standard deviation {^q mn ) and a negative correlation time (— r, which acts as a flag) must be 
specified. These are read in a free format. 

The order of the parameter types in FSN05 is arbitrary; the file is sorted and the time-varying 
stochastic parameters are moved to the top of the file to accommodate the space saving 
implementation of the V p as an nx x nd matrix. 

The stochastic parameters (types 1, 2) are assumed to be zero mean processes. Typically the 
physical process modeled is not zero mean. The non-zero mean is estimated as a constant (types 
4, 5). Thus, the constant (types 4, 5) and the stochastic parameters (types 1, 2) are estimated 
together. The constant can be estimated without a stochastic parameter, but a stochastic 
parameter must be estimated with a constant unless of course the process has a zero mean. 

Since GEODYNH does not have time-varying models the physical processes are modeled by 
estimating constants over consecutive segments of time. For example, the tropospheric scale 
correction in GEODYNII may be modeled over a 24 hour period by estimating a constant over 
the first 12 hours and another constant over the second 12 hours. In OSUORBFS, the one 
constant and a time-varying stochastic parameter would be estimated over the entire span of 24 
hours. This requires the measurement partials from the two consecutive constant estimates in 
GEODYNII to be concatenated. This is controlled by the CONPRT statement. 

OSUORBFS can be implemented as a conventional least-squares sequential estimator by 
specifying all parameters as pseudo epoch state (Type 3) and constant parameters (Types 4, 5). 
The partial file should not be concatenated. The full covariance matrix may be generated. If the 
GEODYNII solution is to be repeated using OSUORBFS in a sequential least-squares step, then 
the full covariance matrix should not be formed. 
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REFSYS 


1 1 1 2 1 3 1 4 1 5 1 6 

REFSYS 910208000000.0000000 

+ 1 + 2 + 3 + 4 + 5 + 6 

COLUMNS FORMAT DESCRIPTION UNITS 


1-6 

A6 

REFSYS - Specifies the True of Reference 
Date (TORD) reference system used by 
GEODYNII. The time MUST be the same as 
the time used on the REFSYS statement in 
FTN05 (G2S). Only TORD is valid. Mean of 
Date (MOD) J2000 is not currently 
implemented 

7 

blank 


21-26 

16 

Year,month,day of reference date (YYMMDD) 

27-30 

14 

Hour, minute of reference date (HHMM) 

31-40 

D10.8 

Seconds of reference date (SS.sssssss) 
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DECORR 


1 1 1 2 1 3 1 4 1 5 1 6 

DECORR 1 0.10 

i 1 1 2 1 3 1 4 1 5 1 6 


COLUMNS FORMAT DESCRIPTION 


UNITS 


1-6 A6 

7 blank 

8 II 


9-10 blank 

11-20 D10.5 


DECORR - controls the computation of the full 
covariance matrix and decorrelation for double 
difference ranges 


= 0, measurements assumed uncorrelated, 
measurements whitened by dividing by the 
GEODYNII supplied weight 
= 1, full covariance matrix computed for 
double differenced ranges, then decorrelated 
and whitened. 


standard deviation for a one-way range meters 
measurement 
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FILSMT 


FILSMT 

1111111111 

-2 + 3 + 4 + 5--- 

O . T. . A . c: 

-- + 6 



Z + .j + ^ t 


COLUMNS FORMAT 

DESCRIPTION 

UNITS 

1-6 

A6 

FILSMT - controls the filter/smoother 
operations, controls the computation of the 
estimates/covariances of the stochastic and 
pseudo epoch parameters (px) and the 
constants (y) 


7 

blank 



8 

11 

= 0, do not filter data 
= 1, filter data 


9 

11 

= 0, do not compute the estimate px at each 
filter step 

= 1, compute the estimate px at each filter step 


10 

11 

= 0, do not compute the covariance px at each 
filter step 

= 1, compute the covariance px at each filter 
step 


11 

11 

= 0, do not compute the estimate y at each 
filter step 

= 1 , compute the estimate y at each filter step 


12 

11 

= 0, do not compute the covariance y at each 
filter step 

= 1, compute the covariance y at each filter 
step 


13 

11 

= 0, do not smooth data 
= 1, smooth data 


14 

11 

= 0, do not compute the estimate px at each 
smoother step 

= 1, compute the estimate px at each smoother 



step 
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= 0, do not compute the covariance px at each 
smoother step 

= 1, compute the covariance px at each 
smoother step 

= 0, do not compute the estimate y at the last 
filter step/first smoother step 
= 1, compute the estimate y at the last filter 
step/first smoother step 

= 0, do not compute the covariance y at the last 
filter step/first smoother step 
= 1, compute the covariance y at the last filter 
step/first smoother step 



UPDTRJ 


+ 1 + 2 + 3 + 4 + 5 + 6 

UPDTRJ 1 

+ 1 + 2 + 3 + 4 + 5 + 6 

COLUMNS FORMAT DESCRIPTION UNITS 


1-6 


7 

8 


A6 UPDTRJ - controls the satellite trajectory 

computation and output in a TORD system to 
file 

blank 

II =0, trajectory is not updated 

= 1, trajectory is updated 
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CONPRT 


+ 1 + --- 

CONPRT 1 

-2 1 3 1 4 1 5 — 

-2 + 3 + 4 + 5-- 

COLUMNS FORMAT 

DESCRIPTION 

1-6 A6 

CONPRT - controls the concatenation of the 
piece-wise measurement partials to the first 
partial location 

7 blank 


8 11 

= 0, no concatenation 
= 1, concatenate 


-- + 6 

-- + 6 

UNITS 
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SATMAS 


+ 1 + 2 + 3 + 4 + 5 + 6 

SATMAS 0 

+ 1 + 2 + 3 + 4 + 5 + 6 


COLUMNS FORMAT DESCRIPTION UNITS 


1-6 A6 SATMAS - controls which satellite number 

system to use 


7 blank 


8 II =0, GEODYN international satellite id 

=1, OSU modified international satellite id 
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c= 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine hsttp (dt , np, nx, npx, nd, ntot , tau, vp, rw, s , rpsm, v, dm, maxnd, 
maxnx , maxnp , maxobs , maxntr , maxntc , pns tdv) 

Apply a sequence of elementary Householder transformations to 
partially triangularize the a posteriori information array, 
this process propagates the SRIF from time t to time t+dt . 

Adapted from ’Factorization Methods for Discrete Sequential 
Estimation' by Gerald J. Bierman, pp. 155-157. 


Version 

9305.1 


Variable 

dt 

np 

nx 

nd 

tau (np) 
Vp ( nx , nd ) 


Comments Pgmr. 

Modified to compute correctly the D. Chadwell 
Rp* matrix. The off-diagonal elements 
were missing from original version. 


Type i/o Description 


r *8 
i*4 

i*4 

i*4 

r*8 

r*8 


Rw ( np ) r * 8 

S (nx+2 *np , ntot ) r *8 


Rpsm(np* (np+1) /2)r*8 


i Propagation interval (secs.) 

i Number of stochastic time-varying 

parameters 

i Number of bias (constant) parameters 

i Number of orbit-related stochastic 

paramters 

i Correlation times 

i The first Nd columns of the Vp matrix 
correspond to the dynamic paramters; 
the last Np-Nd columns are omitted 
because they are in theory zero. 

i Process noise standard deviation 

reciprocals 

i The top Np+Nx rows of S contain the 

SRIF array corresponding to the p 
and x variables; the bottom p rows 
are used to store smoothing related 
terms . 

0 Time-updated array with smoothing- 

related terms stored in the bottom 
portion of S 

1 Upper triangular matrix contains 

smoothing related terms from the 
t-dt to t. 

o Upper triangular matrix contains 

smoothing related terms from t to 
t+dt . 


c 

c 

c 

S on 
Npx 

input : 

S on output : 


c 

c 

c 

c Tp 

Np Nx 

1 Ny I 

1 


I Np | ■'Rp A Rpx 

1 Nx | A Rpx A Rx 

< < 

A Zp I I 1 ~Rp ~Rpx I ~Rpy I 

A Zx I I ==> I ~Rxp ~Rx | ~Rxy 1 

~Zp 1 
~Zx | 

c 

c 

1 Np | 0 0 

i 

1 0 I 

Zw | | I *Rpp *Rpx j *Rpy 1 

*Zp I 

c 

c 

l 


j 


c 

implicit none 
integer*4 

i* j 

, k , 1 , np, nd, nx, npx, ntot , j 1 , j 2 



integer *4 
integer *4 


maxnp , maxobs , maxnt r , maxntc 
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double precision 
double precision 
double precision 
double precision 
double precision 
double precision 
double precision 


z = 0 .dO 
do j=l,np 

if ( tau ( j ) .gt . 0 .dO) then IJRw for 1st order Gauss Markov Model 

dm(j)=dexp( -dt/tau(j) ) 

rw ( j ) =1 .dO/ ( pnstdv (j ) *dsqrt ( 1 .dO-dm ( j ) *dm ( j ) ) ) 
else ! !Rw for random walk model 

dm(j)=l.dO i I If tau. It. 0 flag for random walk 

rw ( j ) =dsqrt { l.dO/dt )/pnstdv(j) 
endi f 
enddo 


vp (maxnx , maxnd ) 

s {maxntr , maxntc) 

v ( maxnp +maxnx+maxobs ) 

rw (maxnp ) , t au (maxnp ) , dm (maxnp ) , z , s igma 

alpha, delta, dt , tmp 

rpsm ( maxnp* (maxnp+1 ) /2 ) 

pnstdv ( maxnp ) 


idiag=0 
idiag2=0 
do j 1 = 1 , np 

idiag2=idiag2+ j 1 
if ( j 1 . le .nd) then 
do i=l,npx 
do k=l,nx 

s (i, 1) =s (i, 1) -s ( i , np+k) *vp (k, jl) 
enddo 
enddo 

if ( j 1 . gt * 1) then 
ict=npx+ j 1-1 
idiag=idiag+ j 1 
index=idiag 
do i=jl,np 
do k = 1 , nx 

rpsm(index) =rpsm(index) -s (ict,np+k) *vp(k, jl) 
enddo 

index=index+i 
enddo 
end if 
endif 

alpha = -rw( jl) *dm( jl) 1! Assumes an uncorrelated process noise 

1 ! COV. 

sigma = alpha*alpha 
do i=l,npx 

v ( i ) =s ( i , 1 ) 

sigma = sigma + v(i)*v{i) 
enddo 

sigma = dsqrt ( sigma ) 
alpha = alpha - sigma 
ict= j 1-1 
index=idiag2 
rpsm (index) = sigma 
sigma = 1 .dO/ (sigma*alpha) 
do j2=2,ntot 
delta=z 

c if ( j2 .eq.ntot)delta=alpha*zw( jl) !! Assume zero mean 

do i=l,npx 

delta = delta + s(i,j2)*v(i) 
enddo 

delta=delta* sigma 
1 = j 2 — 1 


32 



tmp=delta*alpha 
if ( j2 .gt .np) then 
1 = j 2 
else 

iff j 2 . le . (np+1 - j 1 ) ) then 
ict=ict+l 
index= index* ict 
rpsmf index ) = tmp 
end if 
end if 

s (npx* j 1 , 1 ) =tmp 
do i = 1 , npx 

s (i, l)=s(i, j2)+delta*v(i) 
enddo 
enddo 

c s (npx+ j 1 ) , ntot ) =s (npx+ j 1 , ntot ) +delta*zw ( j 1 ) I! Assume zero 

! ! mean 

delta=alpha*rw ( j 1 ) *sigma 
s (npx+ j 1 , j 1 ) =rw ( j 1 ) +delta*alpha 
do i=l,npx 

s ( i , np) =delta^v ( i ) 
enddo 
enddo 
return 
end 
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